
JADAVPUR UNIVERSITY- Course 2019/ 2019 
ELECTRONICS AND TELECOMMUNICATION DEPARTMENT 


Research Report 


Application of Compressed Sensing in Fast MRI Acquisition 


Made by: 
Avrajit Ghosh 



Contents 


1 Research Statement 2 

2 The Compressed Sensing Problem 2 

2.1 LI norm gives the same reconstruction result as LO norm . 3 

2.2 Greedy methods . 3 

3 Application of Compressed Sensing in Magnetic Resonance Imaging 4 

3.1 Analytical transform based Compressed Sensing MRI (CSMRI-LDP) . 5 

3.1.1 Random Sampling reduces aliasing artifacts . 6 

3.1.2 Performance of CSMRI . 6 

3.2 Dictionary learning based Compressed Sensing MRI (DLMRI) . 6 

3.2.1 DLMRI Algorithm . 8 

3.2.2 Performance of DLMRI compared to CSMRI (LDP) . 9 

4 Non-Fourier Random Encoding based CSMRI 10 

5 Future work and Vision in Fast MRI acquisition 10 


1 










1 Research Statement 


Conventional image compression techniques (JPEG) in digital camera first acquires all the samples of 
an image. In the next step, projection of the acquired image in a sparse basis enables it’s compression. 
For example in JPEG-2000, an image of size one megapixel can be represented in as few as 25000 
wavelet coefficients. The decimation in number of coefficients not only allows better storage facility 
for data but operating on fewer number of coefficients consumes less computational time. However 
the question that is to be asked is ” Is it logical for the sensor to acquire so much data only for it to 
be thrown away? Is it possible for the sensor to acquire fewer measurement data ? And if yes can 
the image be reconstructed from these fewer measurements?” In 2004, the answer to this question has 
been addressed by Emmanuel Candes, Justin Romberg and Terence Tao [ ]. They showed that it is 
possible for a signal to be reconstructed with even fewer samples than the sampling theorem given the 
knowledge about signal’s sparsity under certain strict conditions. Today this idea has been applied 
to fields like Magnetic Resonance Imaging, Photography and Network Tomography and is known as 
’Compressed Sensing’. 


2 The Compressed Sensing Problem 

Let x be a sparse signal of dimension N x 1 and A be a matrix of dimension M x N known as the 
measurement matrix. 


y = Ax (1) 

Then signal y is known as the measurement vector which is of dimension M x 1. If x is not sparse, 
then it is possible to project it into a basis where it is sparse say x = Tc. T is an orthonormal basis 
(wavelet, Fourier domain) where the transformed signal (c) is sparse. So the modified equation can 
be written as: 


y = ATc (2) 

It is evident that recovering x or c would require solving an underdetermined system of equations as 
we want to recover N samples from M samples of measurement vector where {M < N ). There are 
infinite possible vectors x that satisfy equation 1 all of which lie on a subspace of rank (N — M). 
However to recover only a specific value of x we utilise the prior assumption that x is sparse. Thus to 
recover x we need to solve the following optimisation problem. 

Minimize) \x\\osubject toy = Ax 

X 

||:c||o is the L0 norm of the vector x which denotes it’s sparsity or the number of non-zero elements of 
the vector. If sparsity of x is S, then to find min ||x||o, we have to check for all (^) possibilities which 
is NP hard to locate the non-zero values in x. So solving the optimisation problem for L0 norm is 
not feasible. Instead of the L0 norm we minimise the LI norm when the measurement matrix follows 
certain conditions [ ]: 

1) Mutual Incoherence 2) Restricted Isometry Property 
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2.1 LI norm gives the same reconstruction result as LO norm 

Due to the NP-hard nature of the optimisation problem, we minimise the LI norm. So the modified 
problem is 


Minimizel Ixl I ^subject toy = Ax 

X 

But how can we guarantee that minimising the LI norm will give us the exact reconstruction results 
that would have been given by minimising the LO norm? Also we ask ourselves the question that will 
the reconstruction result be exact if we minimise the L2 norm. The answer lies in the geometry of the 
Lp balls as shown in Figure-1. 


i^unit: baljs I? ]5 = {x Er R N : Hx|| p < 1} 



Figure 1: Geometry of Lp unit balls for various values of p 

It is evident from Fig-2.ii) that minimising the L2 norm does not give a sparse solution for x. But 
minimising the LI certainly returns a sparse solution for x as shown in Fig-3. The advantage of using 
the LI norm is that the optimisation problem reduces to a convex optimisation problem that can be 
easily solved by the gradient descent method. 

2.2 Greedy methods 

Apart from the LI minimisation problems, there also exists greedy methods for sparse signal recovery. 
Infact these greedy methods provide a faster solution than the LI minimisation problem. Various 
types of greedy methods are: 

1) Matching Pursuit:- The inner product between the measurement vector y is calculated with all 
the columns of measurement matrix A which is denoted by a^. Let a/. =< y,ak > be the coeffi¬ 
cient obtained from the inner product. In the next step, we subtract the effect of that column of A 
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Figure 2: i)Solution subspace of the x (left) ii) Particular solution of x from the subspace that satisfies 
the minimum L2 norm (right). Solution obtained here is not sparse [6] 

which has the maximum inner product. So the updated residue of y in the ith step is represented by 
yi -|_i = yi — OL maXi a rnaXl where a maXi represents the maximum value of inner product obtained at step 
i. This procedure is repeated until we define a maximum number of iteration steps or we define a 
certain threshold for ||yj + i||. The recovered signal is composed of the a maXl ! s where i ranges from 1 
to S, S being the number of iteration steps. 

2) Orthogonal Matching Pursuit:- 

In OMP, the procedure is exactly like MP except the fact that OMP never re-selects the the same 
element a/.. So the updated residue is always orthogonal to all the columns of A i.e a 

3) Compressive Sampling Matching Pursuit (CoSaMP) 

CoSaMP keeps track of an active set of non-zero elements and adds or removes element at each 
iteration. In the first step an S-sparse estimate of x (say xo) is calculated and then the residual error 
||y — Axo|| is projected on the columns of A. The indices of those columns of A which have the k 
largest inner products are selected and are added to the current support of x to get a larger set of 
non-zero elements. An intermediate step is calculated by solving a least square solution and the k 
elements of this intermediate step are used as the new support set of non-zero elements. 

3 Application of Compressed Sensing in Magnetic Resonance Imag¬ 
ing 

Image acquisition speed is an important factor in Magnetic Resonance Imaging as MRI acquisition 
could take a long time and the patient needs to be remain absolutely still while the acquisition is in 
progress. Thus fast imaging technique was in very high demand. As discussed in Section 1 this is a 
classical problem of Compressed Sensing as we want to reconstruct the original image of the brain by 
acquiring fewer samples. 

The data in MRI is acquired in the k-space of the image ( Fourier domain) sequentially over time, 
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Figure 3: Particular solution of x from the subspace that satisfies the minimum LI norm. Solution 
obtained here is sparse. [6] 

However to make the acquisition fast we don’t consider all the samples in the k-space but we un¬ 
dersample the k-space and take fewer k-space measurement data. To recover the image, we exploit 
it’s sparsity in another transform domain where it is sparse. In CSMRI proposed by Lustig et al, a 
definite transform domain is selected and then the theory of Compressed Sensing is applied. CSMRI 
adopts the process of analytical sparsifying transform to make the signal sparse. However it was later 
shown by Saiprasad et al that learning a dictionary for the specific image will provide better sparsity 
and hence improved undersampling rate of the k-space without further degradation of reconstruction 
quality. 

3.1 Analytical transform based Compressed Sensing MRI (CSMRI-LDP) 

In the first step of CSMRI-LDP [ >], the undersampled k-space of the image y is acquired. To recover 
the original image m from the under-sampled k-space F u m we exploit the sparsity of the original image 
by projecting the image onto a sparsifying transform basis if). To guarantee perfect reconstruction by 
the Theory of Compressed Sensing it is required that 

1) the original image have a sparse representation in a transform basis, if) could be wavelet transform 
matrix or finite difference matrix. 

2) Aliasing artifacts will be formed when the k-space of the image is undersampled. The transform 
domain should be chosen such that artifacts are incoherent in that transform domain. 

3) The reconstruction should enforce sparsity in the transform domain and also the reconstructed 
image should be consistent with the original image. 

Accumulating all the three points, the reconstruction problem is then defined as : 
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min \\^m\\isuch that\\F u m — y \\2 <e 

3.1.1 Random Sampling reduces aliasing artifacts 

A simple intuitive explanation has been put forward by Lustig etal which demonstrated the impor¬ 
tance of random sampling in signal reconstruction. An image which is sparse in the spatial domain was 
uniformly undersampled in the k-space. By the property of Fourier Transform, uniform undersampling 
in k-space lead to periodicity in the spatial domain and due to the overlapping of the replicated sparse 
support, the original sparse image could not be reconstructed. 


k-space 


b 


Figure 4: a) Image (wideband) is sparse in spatial domain, k-space of this image is undersampled 
randomly which leads to exact recovery of the the sparse image. 

b) Image is sparse in the wavelet domain. LI norm of the transformed image is minimised in this case. 

[3] 

However when the same k-space of the sparse signal was undersampled randomly, the incoherent 
artifacts due to random undersampling acted as additive noise. It was then possible to recover the 
non-zero values of the image distinctly. 

However in most of the cases the image is not sparse in the spatial domain, so it is projected in known 
basis where it is sparse. Especially when the image is not wideband and is localised in spatial domain 
, the wavelet domain serves as the appropriate basis. So instead of minimising the LI norm of the 
original image m here the LI norm of the transformed image ipm is minimised. 

3.1.2 Performance of CSMRI 

Acceleration of the imaging is proportional to the under-sampling rate of the k-space. It was shown 
that for an 8-fold acceleration the image was exactly recovered for both uniform density and variable 
density sampling. There was a significant improvement in the under-sampling rate of 2.5-3 fold in 
CSMRI. 

3.2 Dictionary learning based Compressed Sensing MRI (DLMRI) 

Compressed Sensing based MRI that is based on analytical transform is limited to an under-sampling 
rate of only 2.5-3. The CS theory states that if the sparsity of signal increases in a transform domain 
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then the signal can be exactly recovered from fewer measurements. To increase the under-sampling 
rate, the transform basis is to be selected in such a way that provides the maximum sparsity for the 
given image. It was shown by Saiprasad and Y.Bresler that adaptive transforms can be sparsify signal 
or images better since they are learnt for the particular image or image class. Moreover exploiting 
patch-based sparsity by dividing the image into several patches can capture local image features more 
effectively . In [ ] a novel CS method has been proposed that learns a dictionary for the image and 
reconstructs the image from highly undersampled k-space data exploiting the sparsity of signal in the 
adaptive transform domain. 

K-SVD dictionary learning [2] is used to learn the dictionary from an image patch (xij)- A sparsity 
constraint T$ is enforced on the coefficients which when linearly combined with the dictionary 
atoms fits the original image patch ( RijX ). So the dictionary learning step can be states as : 


Minimize > ( 
d,t ^ 


| RijX Do/.ij\ 


|o < T 0 




In this equation Rij denotes the image patch of size (y/~N x y/N) extracted from the original image 
of size (N x N). The top- left corner patch of the image is located at pixel (i, j) of the original 
image, = Rijx. t represents the set of of sparse representations which includes all possible 
combinations of vector aij. To enforce the sparsity on transform domain , the number of non-zero 
values of a,j should be less than a specified value of To. To perfectly reconstruct the image patch x 
should satisfy the constraint F u x = y in noiseless case. So in that case where noise is added in the 
k-space the modified optimisation problem is : 


(PO) Minimize y~](| |RjjX - Daij\\ 2 ) 2 + t(IIK - y\\ 2 )s,t\\ a ij\\o < T 0 

D,t l 

ij 

Optimisation problem PO is solved by alternately solving two separate optimisation problems PI and 
P2 till the cost function in PO converges to a certain threshold. Problem PI is the Dictionary learning 
step where K-SVD is utilised to learn the dictionary D using fixed value of x tJ . After the dictionary 
is updated in this step, sparse coding is applied to find the sparse vector a y for all the patches. 

{Pl)Minimize\^(\\RijX - Daij\\ 2 ) 2 s,t\\aij\\ 0 < T 0 \\d k \\2 = 1 
ij 

In Problem P2, the reconstruction is updated by updating x and keeping D, onj fixed. This step takes 
care of the constraint that the k-space of the updated image x is consistent with input under-sampled 
k-space data. 


{P2)Minimize'y^{\\R ij x - Daij\\ 2 ) 2 + 7(11-^ ~ vlU) 

X L ' 
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Oringial MR Image 


k-space of MR Image with added noise 




Figure 5: DLMRI reconstruction with additive k-space noise standard deviation 0.2. Variable den¬ 
sity random sampling mask is used for undersampling of k-space by a factor of 20. [All the re¬ 
sult has been generated from DLMRI software package: http://www.ifp.illinois.edu/ yoram/DLMRI- 
Lab/DLMRI.html ] 

3.2.1 DLMRI Algorithm 

The input to the DLMRI reconstruction algorithm is the undersampled k-space data (y) of the MRI 
image. The algorithm is solved by alternating optimisation problem (PO) AND (PI). The initial value 
of original image x is calculated by taking the IFFT of the zero filled k-space measurement data. The 
next step involves solving the Dictionary Learning Step which finds out the D and a. t j for the initial 
value of x by solving PO. After this step, keeping the dictionary D and sparse representation aiij fixed, 
x is updated by solving (P2). Each pixel value for x is obtained by averaging the contributions of 
patches that cover it. This updated patch averaged value of x is used to recover the data in k-space 
that was initially undersampled. The IFFT of the updated k-space is taken and acts as the input 
of the algorithm in the next iteration. The number of iterations is decided by the stopping criteria 
of the convergence which could be the value of the objective function or the norm of reconstruction 







Norm of difference between reconstruction at sucessive iterations 


45 

T— 

i 40 

1 

1 

X 

■a 

£ 35 

\ 

_ 

* 30 



c 

a> 

0) 

£ 25 



o 

-Q 

a> 

o 20 



0) 

<D 

|t 15 



T3 

a> 

£ io 



*- 

O 

I • 


_ 

z 




0 5 10 15 

Iteration Number 


Figure 6: Norm of difference between reconstructions at successive iterations till convergence 


difference between successive iterations. 

3.2.2 Performance of DLMRI compared to CSMRI (LDP) 

1) With Variable Density Random Sampling (noiseless case): 

In case of Variable random sampling at 5-fold under-sampling, LDP-CSMRI algorithm is unable to 
remove undesirable artifacts. But DLMRI was shown to be free from artifacts and was close to the 
original image. PSNR value for DLMRI after 10 iterations was shown to be 18 dB higher than LDP. 
The CSMRI reconstructed image was shown to have 5% error magnitude in almost half of the image, 
whereas DLMRI reconstructed image was almost free of high magnitude error. 

2) With Cartesian Sampling: 

In the noiseless case, with Cartesian 7.11 fold under-sampling, LDP could not remove aliasing artifacts 
whereas DLMRI acheived a PSNR of 29.5 dB but still had some artifacts. When noise is added in the 
k-space of standard deviation 18.8, PSNR for DLMRI was shown to be 30.68 whereas LDP was still 
unable to remove aliasaing and noise due to zero-filling. 

3) With Pseudo Radial sampling : 

An image of the Lumber spine with 6.09 fold undersampling rate was shown to be recovered by DLMRI 
with less aliasing errors than LDP-CSMRI. 
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PSNR of updated image at each iteration 
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Figure 7: Variation of PSNR of updated image at each iteration step till convergence 


4 Non-Fourier Random Encoding based CSMRI 

The generalised CS MRI reconstruction without any prior sparse basis can be written as : 

min \\ipm\\isuch that\\Em — y \\2 < e 

In the previous sections E was considered to be the undersampled Fourier- Transform matrix( F u ) 
and ip to be the wavelet transform matrix. However encoding the measurement data in Fourier data 
is not well suited for an arbitrary ip as it does not satisfy the RIP conditions for CS reconstructions. 
The choice of random encoding was first proposed by [5] in which it was shown that if E is a random 
Gaussian matrix , then there is a high probability that the RIP will be satisfied for any ip. 

In [ ], three different data acquisition was performed with fixed encoding data size: 1) Random En¬ 
coding, 2) Fourier Encoding 1 (FE1) 3) Fourier Encoding 2 (FE2). It was shown in these experiments 
and results that it has the potential to outperform Fourier-based schemes in high-SNR scenarios. 

5 Future work and Vision in Fast MRI acquisition 

DLMRI has adopted three kinds of sampling mask to demonstrate the effect of under-sampling, Future 
work on MRI compressed sensing can focus on the spatial geometry of the sampling mask that could 
capture most of the information from the k-space of the image. Cartesian sampling mask extracts a 
particular orientation of the original image. One possible solution is to eradicate the effect of a global 
sampling mask by adapting a sampling mask to a certain class of images (MR images). Just like 
the dictionary is adapted according to the original image, the sampling mask can also be adapted to 
extract the frequencies that are predominant in the image class. Using an adaptable sampling mask 
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can result in a higher under-sampling factor than a using a global sampling mask (Cartesian sampling 
mask, Pseudo Radial Sampling mask), 
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